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In the framework of the Fermi-Pasta-Ulam (FPU) model, 
we show a simple method to give an accurate analytical esti- 
mation of the maximal Lyapunov exponent at high energy 
density. The method is based on the computation of the 
mean value of the modulational instability growth rates asso- 
ciated to unstable modes. Moreover, we show that the strong 
stochasticity threshold found in the /3-FPU system is closely 
related to a transition in tangent space: the Lyapunov eigen- 
vector being more localized in space at high energy. 



A large number of theoretical and numerical studies 
have been devoted to the characterization of chaotic high- 
dimensional systems, but in spite of these efforts sev- 
eral fundamental items are not yet fully understood. In 
particular, the relation between Lyapunov analysis and 
other phase space properties like diffusion of orbits, re- 
laxation to equilibrium states, spatial development of in- 
stability remains to be clarified This Rapid Com- 
munication, besides presenting an estimate of the largest 
Lyapunov exponent, is a contribution to clarify this rela- 
tionship in the context of the Fermi-Pasta-Ulam (FPU) 
model 1^. This system was not only the starting point 
of the (re)discovery of the soliton [^, but initiated also 
several studies of the mixed chaotic/ordered phase space 
structure based on resonance overlap criteria B, KAM 
theorem and Nekhoroshev stability estimates M. Statis- 
tical mechanics was also tested on this equation and 
the results showed that ergodicity is not an obvious con- 
sequence of the non existence of analytical first integrals 
of motion. 

There are very few examples of analytical calculations 
of the largest Lyapunov exponent Ai in the chaotic com- 
ponent corresponding to energy equipartition in high- 
dimensional systems In this Rapid Communication, 
we derive an approximate analytical expression for Ai 
which agrees very well with numerical simulations. We 
focus our attention on the asymptotic state of the system, 
when it has reached energy equipartition, i.e. the state 
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where energy is evenly distributed among all Fourier 
modes (relaxation times will not be our concern here). 
One of the main points of this paper is to emphasize the 
relevant role played by some unstable periodic orbits cor- 
responding to Fourier modes. Therefore, we will first de- 
rive the criterion for modulational instability of a plane 
wave on the lattice. 

Denoting by Un{t) the position of the nth atom [n € 
[1, A^]), the equations of motion of the FPU chain read 



(1) 



Mn= Un+1 + Un-1 - 2u„ + 

13 [{Un+1 - Unf^^^ - {Un - V 

where p is an integer greater or equal than 1. We chose 
periodic boundary conditions. Even if the positive pa- 
rameter P can be forgotten by appropriate scaling trans- 
formations of Un , we will keep it in order to make reliable 
comparisons with previous papers, where P = 0.1. For 
sake of simplicity, we consider first the case p — 1 and 
then we generalize to any p-value. 
Looking for plane wave solutions 



(2) 



where 9n(t) = qn — tot and q — 2TTk/N, we obtain 
the dispersion relation uj'^{q) = 4(1 + a) sin^(g/2) where 
a = 12/3(/)q sin^(q/2) takes into account the nonlinearity 
jlOf . The modulational instability of such a plane wave 
is investigated by studying the linearized equation asso- 
ciated to the envelope of the carrier wave. Therefore, one 
introduces an infinitesimal perturbation in the amplitude 
and looks for solutions 

Unit) = [00 + bn{t)] e^^"(*) + [00 + b*it)] e-*^"(*) . (3) 

Introducing this ansatz in Eq. (|l|), after linearization with 
respect to 6„ but keeping the second derivative (contrary 
to what has been done for Klein-Gordon type equation 
[]lT|), we obtain 

bn - 2iujbn^ (1 + 2a)(6„+ie*« + 6„_ie-*'? - 2cos(g) 6„) 
+ K^i - 2 cos(g) bn) (4) 

Assuming further 6„ = A e^W"""*) + B e-*(Q"-"*), we 
finally obtain the following dispersion relation 



4(1 + 2a) sin^ 



(f] - Lo)2 - 4(1 + 2a) sin^ ' ^ 




4a^ (cosQ — cosq)^ (5) 
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This equation has 4 different solutions once q (wavevec- 
tor of the unperturbed wave) and Q (wavevector of the 
perturbation) are given. If one of the solutions is complex 
we have an instability of one of the modes (g± Q) with a 
growth rate equal to the imaginary part of the solution. 
Therefore, one can compute the instability threshold for 
any initial linear wave, i.e. any wavevector and any am- 
plitude. For example for q — Q we find that the solution 
is obviously stable since the zero-mode, corresponding to 
translation invariance, is completely decoupled from the 
others. For q = n, the expression for the growth rate is 

r(7r, Q)= 2(^(1 + a)(4 + 8a) cos2(g/2) + cos4(g/2) 



-l-a-(l + 2a) cos^(g/2) 



1/2 



(6) 



A simple analysis of this function shows that the first 
unstable mode is the nearest mode corresponding to Q = 
2tt/N. Computing the critical value of the parameter a 
above which r(7r, 2'k/N) is positive, we obtain the critical 
energy Ec for the 7r-mode. It reads 



Er 



2iV . 2.7r 7cos2(f )-l 
9/3 ''''' ^N'Scos^if)-! 



(7) 



This analytical expression is in agreement with the previ- 
ous p^-p^ approximate expression Ec — Tr^/SiV/? valid 
only in the large N limit. Above this energy thresh- 
old, the vr-mode is therefore unstable and gives rise to a 
chaotic localized breather-like excitation, able to move 
very fast in the system, collecting energy from high- 
wavevectors phonons until its disappearance leads to the 
final energy equipartition [l5|] . 

When the energy increases, the region of instability 
extends to a larger region of wavevectors, and in par- 
ticular the most unstable mode Qmax increases until the 

asymptotic value Qmax — 2 SiTCCosi^^J 8/\/3 — 4) ~ 0.427r 
is reached. It is important to notice that, for sufficiently 
high energy, the rescaled growth rate t/t{'k, Qmax) does 
not depend on the energy density. The growth rate is 
plotted in Fig. |^ at high energy. 

We performed some simulations of the system with a 
6th-order symplectic integration scheme adopting as ini- 
tial condition the 7r-mode and computing the Lyapunov 
exponents after the transition to equipartition. We used 
the algorithm proposed by Benettin et al , where the 
full set of tangent vectors is periodically reorthonormal- 
ized using the Gram-Schmidt method; the Lyapunov ex- 
ponents are then obtained from the time average of the 
logarithms of the normalization factors. The results are 
plotted in Fig. ||. We have also checked our results by 
performing the numerical integration directly in Fourier 
space, paying particular attention to the Lyapunov eigen- 
vectors. 

Let us present now the analytical estimation of the 
maximal Lyapunov. As the system is symplectic, the 



usual pairing rule is valid and moreover Pesin's theorem 
allows us to identify the Kolmogorov-Sinai entropy of the 
system with the sum of all positive Lyapunov exponents. 
As the spectrum was shown [ |l7t to be approximately 
linear at high energy (see the inset of Fig. ^), one can 
relate the Kolmogorov-Sinai entropy Sks with the max- 
imal Lyapunov exponent, namely 



N 



Let us define the instability entropy 

N/2 



(8) 



(9) 



where the sum is over all positive growth rates | fi8[ . 
The crucial physical hypothesis of this paper is that 
Sks — Sie{t^), we then obtain the following analytical 
expression for the maximal Lyapunov exponent: 



N/2 



(10) 



Using the expression (|^), we can then compute the max- 
imal Lyapunov exponent. Fig. |^ attests that the analyt- 
ical expression (^^ is very accurate. In the same figure 
the data obtained with a completely different approach, 
developed by Casetti, Livi and Pcttini (CLP) 1^, are also 
shown. The two methods give almost identical results, 
apart at very low energy, where the CLP-findings is 
in better agreement with our numerical data. However, 
our approach is definetely simpler and relies on the anal- 
ysis of unstable periodic orbits, while the CLP-one on 
Riemannian differential geometry. 



It is remarkable to note that Chirikov |19| found sim- 
ilarly the maximal Lyapunov exponent of the standard 
map at high energy by averaging over the phase space 
the maximal eigenvalue associated to the main hyper- 
bolic point. It corresponds in our case to averaging the 
growth rate (^ for the unstable periodic orbit q — n over 
the equilibrium equipartition state (where all modes have 
the same weight). A similar approach is known as Toda 
criterion |Q and although it cannot be used as a signa- 
ture of chaos, it can give an approximate estimation of 
Al. 

In fact, one can understand this average in a better 
way by recalling that the modes {7r/2}, {27r/3}, {tt} cor- 
respond to the simplest unstable periodic orbits and are 
also the only three one-mode solutions of the /J-FPU 
problem. The calculation of the instability entropies of 
this three modes shows that they are extremely close one 
to another contrary to the value for other modes. A cor- 
rect approach would be to apply the zeta-function for- 
malism pl[ to this system, if feasible. 
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At high energy, expression (g) can be simphfied as 
T(7r, Q) ~ ^/a f{Q) where f{Q) is energy independent. 
Therefore the growth rate scales with the amplitude (j)o, 
and as E = N (80q + 64/3(^q), it means that the growth 
rate and therefore the maximal Lyapunov Ai scales with 
(E/Ny/'^ at high energy. This result is in contradiction 
with Ref. ||2^] but in agreement with Ref. 

In fact, a similar approach gives also very good re- 
sults for other powers 2{p + 1) in the coupling potential. 
The expression of the growth rate is then the same if we 
use a = /?igtlJ(2<^o sin(g/2))'^ Fig. | shows that the 
results are once more in very good agreement with nu- 
merical estimates. One derives easily that the maximal 

1 1 

Lyapunov scales at high energy like Ai ^ {E/N) ^ 2(p+i) 
psf . It is important to stress that in the limit of hard po- 
tential {p — > oo) we find the exponent 1/2, analogously 
to billiards [ p4| . In the low energy limit, however, all 
models have the same scaling behaviors as expected. 

Plotting \i{E/N) in a log-log scale we observe (see 
Fig. 1^) that the two asymptotic linear behaviors are sep- 
arated by a knee at intermediate energy density. An es- 
timation of this transition region can be derived assum- 
ing that the linear and nonlinear contributions to uj^q) 
should be of the same order. We obtain a ~ 1 i.e., an 
energy density of the order of 1//3 (this is equivalent to 
the estimation given by mode overlap criterion Q). 

This knee corresponds to a stochasticity threshold 0] 
which defines the crossing from weak to strong chaos. 
But we have also found that it corresponds to an inter- 
esting transition in tangent space. Considering the nor- 
malized Lyapunov vector Vi associated to the maximal 
Lyapunov Ai, we can introduce the participation ratio 
PI 



/ N 




(11) 



where the first (resp. last) N components of Vi are asso- 
ciated to the evolution of linear perturbation of w„ (resp. 
iin) in tangent space. The quantity (11) has been used 
in different contexts and for example in dynamical sys- 
tems as an indicator of localization: it is of order N 
if the vector is extended and of order one if localized. We 
have found that the stochasticity threshold corresponds 
to a crossover from an extended state in tangent space 
to a more localized state, as attested by Fig. Q. The two 
examples of Lyapunov vectors reported in Fig. ^ show 
the clear difference between low and high energy den- 
sity. This transition is very reminiscent of the metallic- 
insulating transition in finite samples |27|. 

We conclude by stressing again that we have computed 
the maximal Lyapunov for a high dimensional Hamilto- 
nian system by using a simple analytical approach, based 
upon modulation instability analysis of linear waves. The 
results obtained here are in excellent agreement with our 



computer simulation results. The success of this cal- 
culation suggests that this Lyapunov estimation could 
be extended to other high-dimensional Hamiltonian sys- 
tems. Moreover, we have shown that the strong stochas- 
ticity threshold Q is not a threshold to energy equipar- 
tition, since equipartition can be always obtained, al- 
though on longer time scales It rather corresponds 
to a crossover from extended to more localized state in 
tangent space. 
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FIG. 1. Shape of the growth rate r(7r, Q)/r(7r, Qmax) for 
sufficiently high energy. The diamonds corresponds to the 
case A'' — 256 and the solid curve to the asymptotic shape 
when A'^ goes to infinity. 
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FIG. 2. Comparison of the analytical estimate with nu- 
merical results for the maximal Lyapunov exponent. The 
solid curve corresponds to our estimation (Eq. llOl) , the dashed 
curve to the estimate using Riemannian differential geometry 
(see Ref. [^) and the triangles to our numerical results (for 
the /3-FPU, i.e. p = 1). The dotted curve (resp. diamonds) 
corresponds to the analytical estimate (resp. numerical re- 
sults) in the case of a power p — 2. In the inset, we plot the 
N positive Lyapunov in the case E/N — 4200 and p = 1. 
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FIG. 3. Participation ratio ^ versus the density of En- 
ergy E/N, for different lattice sizes; diamonds for N = 64, 
triangles for N — 256 and squares for N — 1024. 
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FIG. 4. Localization in tangent space of the Lyapunov 
vector for N — 256. The dashed curve corresponds to a 
generic localized Lyapunov vector at high energy density 
E/N ~ 2.5 10^, whereas the solid curve (shifted by -0.2) cor- 
responds to a generic delocalized one at low energy density 
E/N ~ 2. 
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